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ABSTRACT 


Changes in the globally integrated absolute angular momentum of the atmosphere 
were computed from the Fleet Numerical Oceanography Center NOGAPS wind analyses 
and compared to astronomically measured changes in length of dav (LOD) obtained 
from the U.S. Naval Observatory, Washington D.C. The two time series were subjected 
to both time and frequency domain analysis. In the tine domain, digital filters were used 
to isolate seasonal and subseasonal components. In the frequency domain, energy den- 
sity, coherence and phase were computed over periods from 2 days to 1000 days. Over 
90% of the total variance in astronomically determined LOD can be explained by 
ineteorological phenomena. Fluctuations in LOD are coherent and in phase with. fluc- 
tuations in the globally integrated angular momentum of the Earth’s shell (crust, mantle 
and oceans; liquid core is excluded) at almost all periods less than 365 days. Annual 
fluctuations in LOD appear to originate in the midlatitudes and propagate equatorward. 


Subseasonal fluctuations (30 to 100 day periods) appear to be a tropical phenomena. 
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I. INTRODUCTION AND THEORY 


A. BACKGROUND 

Fluctuations in the Earth’s length of day (LOD) were surmised by carly astronomers 
who noted cumulative irregularities in the Earth’s expected latitudinal position relative 
to daily star transits. Newtonian physics explained only the tidal component of these 
short period changes in the Earth's angular momentum. In the late 1930's, pendulum 
clocks were accurate enough to measure daily star transits to within | ms.! Stoyko (1937, 
cited in Munk and MacDonald, 1960) reported a 2 ms cumulative seasonal variation in 
LOD. He noted that the length of day measured in July was longer than that measured 
in May by 2 ms. Atomic clocks, developed in the 1950's, improved time keeping accu- 
racy by several orders of magnitude. Present day cesium clocks have accuracies better 
than 10-* s per day (or .00001 ms) (Ramsey, 1988). At the same time, improvements in 
astronomical measurement methods led to more precise determination of star and quasar 
transits. 

A number of studies have been undertaken to link the observed changes in LOD to 
meteorological phenomena. Starr (1948, cited in Munk and MacDonald, 1960) pro- 
posed that the observed changes in the Earth’s LOD occurred because the solid Earth 
and atmosphere exchanged angular momentum on a seasonal basis. Madden and Julian 
(1971) found 40-50 day oscillations in pressure and zonal winds centered in the tropics, 
and Lambeck and Cazenave (1974) subsequently linked these atmospheric oscillations 
to changes in LOD. Lambeck (1980) identified the amplitudes and phases of the seasonal 
changes and listed other, non meteorological components which could be responsible for 
LOD changes. Rosen and Salstein (1983) correlated short period (less than seasonal) 
henuspherical changes of atmospheric angular momentum and LOD with polar and 
subtropical jet streams. Eubanks et al. (1985) performed a spectral analysis of atmo- 
spherically derived angular momentum with a composite of astronomic, satellite and 
Lunar Laser Ranging (LLR)2 LOD data. Barnes er al. (1983) extensively studied 


1 The best gravity pendulum clocks have accuracies within | ms per day or one part in 108 
(Ramsey, 1988). This is barely accurate enough to detect the fairly large seasonal changes in the 
LOD. 


2 The astronauts left a corner reflector on the moon which reflects a LASER signal back for 
extremely accurate short term measurement of Earth rotation parameters. 


changes in LOD over the period 1979-1981. Morgan er al. (1985) used least squares 
analysis to study the seasonal as well as subseasonal nature of oscillations in the LOD 
from 1981 to 1984. Swinbank (1985) examined the torques which produce changes in 
LOD. Wolf and Snuth (1987), and Eubanks er al. (1986) compared changes in atmo- 
spheric angular momentum with astronomically determined LOD during the El Nino of 
1983. 

Most studies have examined atmospheric angular momentum and LOD using time 
mieasurements from the Bureau International de lITleure (BIH), augmented by LLR and 
other precise rotation measurement methods, to compare with gridded zonal wind data 
from the U.S. National Meteorological Center (NMC) or the European Centre for Me- 
dium Range Weather Forecasts (ECMWF). This paper analyzes the wind-derived at- 
mospheric angular momentum budget for the period January 1983 to July 1987 and the 
simultaneous changes in the Earth’s LOD using time data from the U.S. Naval Observ- 
atory (USNO) and zonal wind data from the U.S. Navy Fleet Numerical Oceanography 
Center NOGAPS model.3 The atmospheric angular momentum changes are then exam- 
ined over two latitude bands to investigate the horizontal structure of the seasonal and 
subseasonal components of the atmospheric angular momentum and contributions to 
theme oe 


B. ASTRONOMIC LENGTH OF DAY 

Astronomic length of day is measured by comparing elapsed time for meridional 
transits of stars or quasars with a reference time determined by atomic clocks. The time 
it takes for exactly one rotation of the Earth (360 degrees of longitude change) with re- 
spect to the fixed, inertial reference system of distant stars is a sidereal day with a mean 
value of 86,164.0905 seconds. Earth time is generally kept with respect to rotation about 
the sun: the solar day. A mean solar day is slightly longer and subject to more variability 
than the sidereal day because the Earth is in motion about the sun. Time for a solar day 
(LOD,) is defined as exactly 86,400 seconds. LOD is defined as the actual elapsed time 
(measured by cesium clocks) for a meridional transit of a point on Earth with respect to 
an inertial frame. Since this is with respect to sidereal time, a conversion 1s applied to 


make it consistent with solar time keeping. 


3 NOGAPS is Naval Operational Global Atmospheric Prediction System. It is a mathematical 
model which generates a global wind and pressure analysis on a 2 % by 2 “% degree gnd spacing. 


LOD 1s related to the Earth rotation rate (w) by: 


LOD, d(UT1) 


ee ae an) UW) 





where LOD and UT! are both measures of the actual elapsed time of the solar day, 2 
is the mean sidereal rotation rate (7.2921 x 10-° s-') of the Earth, and /AT is the precise 
time reference determined from atomic clocks. UT1 is an angular measurement of 
Earth’s polar rotation with respect to sidereal time passage which 1s usually converted 
to tine. 


Changes in the Earth’s rotation rate with respect to accepted standard values are:4 


LOD LOD a(UT\ — IAT. 
0M = Q > = © _ = ae! aie) 


LOD LO Dea dt : (2) 





Taking the difference between LOD and the standard, LOD,,, to be 0OLOD produces: 


OO ape LO Dips GET! = 1 Al aa) 


ee ee oy di : (3) 


where LOD,, is taken equal to LOD, and JAT,,, 1s taken equal to a sidereal day, although 
other conventions could be used with the relationship sull holding (Eubanks er al., 1985). 
The approximation in (3) 1s valid because, using a time step of one sidereal dav, LOD,, 
the ratio of dbLOD to LOD 1s the same as the ratio of bLOD to LOD, to within 10-'s. 
Using finite differences to approximate the UT1 derivative: 


6LOD = - Lo,{ (4) 


(i 1A Tyea)e Sa ciel eid Laat ; 
At 


Eq. 4 relates daily observations of star and quasar transits directly with the precise 
elapsed time determined by atomic clocks. I[t 1s computed on a daily basis (Ar = | day) 


from astronomic data gathered all over the world to help eliminate single site bias. This 


4 By subtracting standard values, which are the international definitions of the solar and 
sidereal days, the remaining quantity represents a change from the standard as well as a change with 
respect to exactly one solar day. The actual measured lJength of day is taken with reference to a 
beginning epoch (1 January 1900, by convention). Thus the actual length of day (LOD) 1s the 
running sum of the changes in LOD from the beginning of the epoch. Thus, the term dLOD, as 
developed in Eq. 4 represents a tume derivative only in the sense of it being the first difference of 
the LOD time sens. If the first difference were not used, the numbers associated with the Iength 
of day would be so large as to be cumbersome. 


quantity 1s Compared with the atmospherically derived equivalent, 6LOD,,.., developed 


atm? 


in the next section. 


C. ATMOSPHERIC ANGULAR MOMENTUM AND LENGTH OF DAY 

Over long periods of time (centuries and longer) Earth angular momentum (Jf,,,n) 
slowly transfers to the moon through solid body tidal friction; hence LOD increases 
slowly and Mj,., decreases. For periods less than 10 years, this effect is negligible; \/-,4, 
is essentially conserved (Eubanks ez al.; 1985, Lambeck, 1980). Thus, changes in M,,.,, 


denoted by O12 7. villsunmte zero: 


OMearih = OM jen + OM jcean + OMe. + 6M = Q. (3) 


core 


Since 06M.,,. and OM,,., Interact in an oscillatory manner with a period of 10 years or 
greater, 0.\/.,,, appears in short data records (under five years) as a linear trend and is 
easily removed by detrending. The three remaining momentum elements: 60/,,,,, OV nen 
and 6M,..., interact to produce changes in the Earth measured LOD. 6M,..,, 18 difficult 
to measure directly. To exchange momentum with the rotating Earth, the oceans must 
have unimpeded zonal flow about the globe. Only the Antarctic Circumpolar Current 
(ACC) has such characteristics. Lambeck (1980) agrees with the earlier findings of Munk 
and MacDonald (1960) that, at most, the ACC could have a semiannual amplitude 5% 
of the atmospheric amplitudes. Thus, with 6M,,,, = 0, and disregarding inputs from 


the ocean and core: 


OMsnen = — OMam- (6) 


The atmosphere of the Earth represents about 10-° of the total mass of the solid 
Earth (solid Earth includes crust, mantle, oceans and atmosphere; everything except the 
liquid core). Because of the large difference in atmospheric and Earth masses, relatively 
large changes in the zonal velocity field of the atmosphere are required to produce small 
Changesmmm yen 


The angular velocity (w) of the Earth’s shell (solid Earth minus the atmosphere) 1s: 


M shett 


a (7) 


l shell 


where M,,., is the angular momentum of the Earth’s shell and /,,., 1s the polar moment 


of inertia. Again, after subtracting standard values from 6M,,.., changes in w are: 


; OM 
iia shell . 


(8) 
l helt 


Overbars indicate average values with respect to time. Since /,,., changes little on time 
scales less than decades, it is taken as a constant: 7.103 x 10° kgm. Substituting 


— 6M,,, from (6) into (8) and combining with (3) yields: 


LOD. OMe 
ey (9) 
Herel 





The assumptions leading to Eq. (9) imply that fluctuations in the length of day are di- 
rectly proportional to fluctuations in the atmospheric angular momentum. One purpose 
of this thesis is to see if independent observations of dLOD and 6 M.,,,, are indeed related 
as predicted by (9). 

The atmospheric angular inomentum at any time can be measured by evaluating a 
volume integral of the east-west (wu) component of the wind from the surface of the Earth 
to the top of the atmosphere. 


2% ae r=00 
We. = \, [ofc pU(r, 6, 6) r° cos) dr dé dO, (10) 
Z 


where r is distance from the center of the Earth and U is the absolute zonal velocity, 
U = rQcos@+u, and u the zonal wind. The hydrostatic relationship allows pressure 
coordinates to be substituted for p dr: 


cl 


dr = — Zp (11) 


Since the depth of the atmosphere is very small compared to the radius of the Earth, r 
is accurately approximated by the mean radius of the Earth, R, (6.37 x 106m). With the 
above substitutions into (10), M,,,, becomes: 


Imb 


27 — sfc 
More = a | | (REQ cos(b) + u)R3 cos*() dp dob dO. (12) 
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Separating the two parts of the integral and simplifying, 


M 


atm — g 





2nRe p1000mb p> ; 
| oil [is cos*($)[u] db dp 
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Re p2n 3 
+o [| _ cos'(b) AP do dd. (13) 
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The first term is the wind or motion term (Barnes et al, 1983). The zonal wind av- 
eraged by latitude bands is indicated by [u] . Integration limits of 100 mb and 1000 mb 
for the motion term were selected because-the NOGAPS wind analysis extends to only 
100 mb and because use of surface winds would add unnecessary complexity to the at- 
mospheric angular momentum calculations. Since winds above 100 mb constitute less 
than 10% of the mass of the atmosphere, exclusion of the stratospheric winds should 
introduce fairly small errors. Rosen and Salstein (1985) found the stratospheric contrib- 
ution to VW, to be positive at semiannual periods and negative at annual periods. In 
other words, the stratosphere opposes the troposphere (180° out of phase) at annual 
cycles, and adds to tropospheric angular momentum at semiannual cycles. The ampli- 
tudes of these signals are less than 10% of the tropospheric term for annual periods, and 
30% of the tropospheric term for the semiannual term. Surface winds are not used be- 
cause of conflicts where surface pressure is near 1000 mb (most of the globe). In that 
case the near surface winds would be doubly represented. 

The second term of (13) is the pressure, or matter, term (Barnes et al., 1983). Most 
studies to date have dropped the pressure term for lack of accurate surface pressure data. 
Barnes er al. (1983) used the special observations from the GARP (Global Atmospheric 
Research Program) of 1979 to evaluate the pressure term and its contribution to M,._. 
They concluded that pressure variations in M,,,, are in phase with the wind term and 


fairly constant at 10% of the wind term amplitude. Indirect studies by Morgan ev al. 


(1985) show the pressure term to be 10% of the wind term at seasonal periods and no 
more than 25% at shorter periods. 

Standard values are subtracted from .\/,,,, to give OAZ,,,,. The standard value for both 
the wind and pressure terms are taken to be Zero. This is consistent with the atmosphere 
having a net angular momentum of zero with respect to the Earth. Neglecting the 
pressure term in (13) and using (9), changes in the Earth observed angular momentum 


caused bv the atmosphere are: 


T 
i 





LOD pp aa. 1h) db d 14 
$LODam = SE ia [fxd cos (b) dep dp. (14) 


Eq. 14 represents the change, in units of time, that changes in the zonal wind will have 
on the rotation of the Earth. Thus, the independently determined measures of 0LOD 
from (4) and (14) should be the same if the atmosphere ts, in fact, responsible for forcing 
short period changes in Earth observed length of day. The data processing and the 
variance computations are described in the next chapter, and all of the results, impli- 


cations and interpretations are given in Chapter III. 


Il. DATA ACQUISITION AND PROCESSING 


A. U.S. NAVAL OBSERVATORY LOD MEASUREMENTS 

Measures of OLOD for the period January 1983 to July 1987 were obtained from the 
U.S. Naval Observatory, Washington, D.C. (USNQO). Included in the data set were both 
the Naval Observatory and the French Bureau International de l’Heure (BIH) estimates 
of bLOD. For notation purposes, 0LODysxo will refer to the USNO data and 6LOD,,, 
to the Frencheq@aca 


B. NOGAPS WIND DATA 

Twice daily analyses of the global zonal winds for the January 1983 through August 
1987 observation period were obtained from the NOGAPS model at the U.S. Navy Fleet 
Numerical Oceanography Center, Monterey CA. The data were in the Navy Environ- 
mental Data Network (NEDN) format. Reporting levels for the winds were standard 
levels plus 925 and 250mb.5 The 2 % ° by 2 % ° grid produces 10,512 reporting points 
for the winds. A FORTRAN program was written to decode the NEDN format, ma- 
nipulate the 11 by 10,512 element matrix and reduce the data to a single value of 
OLOD sm. Eq. 14 was approximated by zonally averaging the matrix, integrating verti- 
cally (using a trapezoidal approximation for the unevenly spaced pressure levels), ap- 
plying the cos? latitude terms and summing over latitude. Since the integration over 
latitude was performed last, the program was easily modified to get d6LOD,,,, values for 
specific latitude bands as well as the globe as a whole. Appendix A contains the source 
code for the program. 

The archived data had few (less than 4%) missing data points. All of the missing 
observations in the wind data were interspersed such that linear extrapolation of the 
mussing Winds could be made. The resulting time series was smoothed with a three-day 
moving average filter and decimated to 1622 elements. 

Similar steps were taken in determining OLOD,,,, values for latitude bands within 
15° and 20° of the equator, respectively. Values for the remainder of the globe were ob- 


tained by subtracting the latitude band series from the original 6LOD,,,, series. The 


5 Standard levels include surface, 1000, 850, 700, 500, 400, 300, 200, 150 and 100 mb. Surface 
winds were not used to calculate 6LOD,,,, because of the possible confusion with the 1000 mb 
winds. 


purpose of the regional analysis was to investigate the possible physical mechanisms 


contributing to the LOD and ,,,, relationships. 


C. ANALYSIS TECHNIQUES 

Fig. | on page 22 shows the global 6LOD,,,, and the 0LOD,.., plotted together. The 
time data have an obvious decreasing linear trend, while the atmospheric data do not. 
The means of both series were removed prior to plotting. The trend in the time data ts 
from decade fluctuations in the LOD caused by core/shell coupling (Munk and 
MacDonald, 1960 and Lambeck, 1980). The trend was removed by subjecting both series 
to a least squares algorithm found in Bendat and Piersol (1971, equation 9.20). The two 
detrended series are plotted in Fig. 2 on page 22. All of the series developed in the 
analvsis were detrended with this scheme. 

1. Time Domain Analysis 

Both series were subjected to linear filtering to isolate seasonal and subseasonal 

components. First, the high frequency noise and tidal oscillations were removed with a 
sine squared weighted moving average filter (also known as a Hanning filter) of 23 points 
(Fig. 3 on page 23). Weights for such a filter are determined by Robinson and Silvia, 
fie? S): 


wh 


WV, = a) 


n= 1,2,3..., M, (15) 
Where }V, are the filter weights and M 1s the total number of weights. The sin? weights 
minimize spectral distortion through unwanted side lobe leakage at the expense of a 
slightly smeared frequency spectrum (Robinson and Silvia, 1978). By keeping the filter 
symmetric, phase information of both series is preserved. Symmetric application of the 
filter results in the general expression of the transformed series y,, (Robinson and Silvia, 
mae O14): 
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The frequency response, or convolution, of a symmetric, moving average filter of M 


points 1s given by: 


An 


| 
HY) = — (hg + 2) hy cos(2n xf), (17) 
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and fh, is the center weight of the full series )f. The term A, will always equal one, since 


) weights of the symmetric series HV, 


it is the center weight. Subtraction moving average filters have a transfer function of 
| — Hi) and serve to pass high frequencies. Sine squared weights produce transfer 
functions almost devoid of side lobes. Plots of H(/) for each of the filters used are con- 
tained in Appendix C. 

One drawback to using convolution techniques is the loss of Af — 1 data points 
(half at each end of the series) with each pass of the filter. This is particularily noticeable 
when using a long string of filter weights often required for lowpass systems. Nonethe- 
less, the record length of the data (1622 points), affords sufficient length to use fairly 
long series of weighted moving average filters. 

To isolate the seasonal terms in the 0LOD,..,5 and dLOD,,,, series, a 75-point 
sin? weighted lowpass filter was applied twice to both series. The effect of successive fil- 
terings in the time domain 1s multiplication of the transfer functions in the frequency 
domain. The transition band with a double pass is narrower than that of a single pass. 
Periods less than 50 days are completely attenuated. The amplitude of annual periods 
are attenuated by 5% and the semiannual terms are diminished by 18%. The resulting 
series are plotted in Fig. 4 on page 23 and the transfer function is plotted in Fig. 18 on 
page 47. To further separate the seasonal terms into annual and semiannual compo- 
nents, both series were subjected to a 183-point unweighted moving average filter. The 
sin? weights would have required too many points to adequately filter out the semiannual 
terms. Using the “boxcar” unweighted filter of 183 points put the first zero at exactly the 
semuannual period, allowing almost complete attenuation of those signals. Disadvan- 
tages of this type of filter are the fairly large side lobes at frequencies above the first zero 
crossing (183 days). Since the boxcar has a fairly narrow transition band, the remaining 
annual terms were still about 65% of their original amplitude. The filter response is 
contained in Fig. 19 on page 48 in Appendix C, and the filtered series is Fig. 5 on page 
24. To isolate the semiannual terms, the series obtained from the 183-point boxcar filter 


was subtracted from the twice run 75-point weighted series. The resulting semiannual 
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series Was attenuated less than 15% at periods of 200 days. Response for this series 1s 
contained in Fig. 20 on page 48. 

Oscillations with periods from 20 to 100 (subseasonal) days were isolated with 
high and lowpass filters. The highpass filter consisted of two passes of a subtracted 
moving average of 75 points. Lowpass was achieved with a single run of a 23-point 
weighted moving average filter. The resulting series are plotted in Fig. 7 on page 25, 
Fig. 21 in Appendix C contains the frequency response of the filter. 

For latitude band comparisons, the 0LODy.x9 Series Was separated into rela- 
tively high frequency (subseasonal) and relatively low frequency (seasonal) components. 
The transition band was from 1/60 day-' to 1/150 day-'. These two series (heavy lines 
Mees eceand /), OO Dy concen 290 OLOD -owocunseass Were then individually compared 
to the corresponding OLOD,,,, time series generated from tropical and extratropical zonal 
winds. These plots are contained in Figs. 8 through 11. Namung conventions for these 


series are: 


© dLOD ysxo-seasonan. Seasonal (150 to 365 day periods) values of USNO determined 
Pep 


© dLODgsxo-surseas. DUOSeasonal (20 to 100 day periods) values of USNO determined 
OLOD. 


© dLOD,._\5: Atmospherically derived 6LOD from 15°N to 15°S. 


®© §LOD wes: Atmospherically derived 6LOD from 15°N to the North Pole and 15°S 
fomrie South Pole. 


© dLOD m2 ANd OLOD ym: Atmospherically derived dLOD with cutoff at 
20°N and S. 


2. Frequency Domain Analysis 

Prior to. Fourier analysis, each series was tapered with a cosine bell to reduce 
leakage at the ends of the time series data. Both series were extended to a length of 2048 
points by the addition of zeros to each end of the record. The effects of the tapering were 
removed after the transform by dividing each element by the ratio of the area of the ta- 
per (the cosine) to the rectangular boxcar area replaced by the taper; a factor of 0.875 
(Bendat and Piersol, 1971). | 

A record length of 2048 permitted the use of the Fast Fourier Transform (FFT) 
algorithm (Bendat and Piersol, 1971), significantly cutting down the computation time 


required for the transforms. After transforming, the complex valued cross spectral am- 


plitude coefficients were determined from the complex amplitude coefficients of the time 


and atmospheric series (X,,, and Y,,,). Cross spectral coefficients are: 


mS Boo. =e (18) 


where » represents the index number of the +. +1 resolvable frequencies from a sample 
size of V and X* ts the complex conjugate of X. The three complex valued sens 
X, Y, and, XY, were then smoothed in order to increase the confidence intervals of the 
power spectral density functions as well as the coherence squared (y?) and phase (@) 
functions. Prior to smoothing, the spectral estimates were prewhitened by multiplication 
by the square of the frequency at each estimate. Jenkins and Watt (1968) state that es- 
timates of coherency and phase should be based on white spectra. This technique (also 
used by Eubanks er al., 1985), compensates for the f? power law followed by both the 
time and atmospheric data. The complex spectral and co-spectral estimates were 
smoothed with a seven-point unweighted moving average to give 14 degrees of freedom 
for calculating confidence and significance intervals. Energy density estimates were 


formed by taking the magnitude of the complex coefficients: 


oe a 


Coherence squared is defined as: 


‘ 2 
IXiny* Yiny| 


| (20) 
G{n)G,(n) 
Phase is determined by: 
[rey 2.5] 
¢ = ta’ (21) 
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Eqs. 19, 20 and 21 were calculated and plotted for the 6LODysno and OLOD,,, series, as 
well as the tropical and extratropical 6LOD,,, series against the ODLODysno . Figs. 12 


through 16 are the plots for these series. 


3. Harmonic Analysis 
For comparison with earlier findings of Rosen and Salstein (1985), Morgan es 
al, (1985), and Langley er a/. (1981), a harmonic analvsis was performed to obtain the 
annual and semiannual terms of both time and atmospheric signals. The Fourier anal- 
ysis was done for the three calendar years for which the data were continuous: 1984, 


[985 and 1986. Phase and amplitude were determined from the sine and cosine sums of: 


? 365 ; > 365 


Where f, is the frequency of interest (1/365 day" and 1/182.5 day"') and ¢ is the integer 
time in days from the beginning to the end of the year. Amplitude for the frequencies is 
J Adios + Abo and phase is tan-'(A,,;,/A,..,) converted to davs from | January of that year. 
Results and comparisons with earlier findings are in Table 3. 


Il. RESULTS 


A. GLOBAL COMPARISON IN THE TIME DOMAIN 
The detrended 0LOD, x. and O0LOD,,, series (Figs. 2 and 3) match fairly closely. 
The linear correlation coefficient for the two series with tides and high frequency noise 
removed is 0.990.6 From Fig. 2 on page 22 it can be seen that the signals with the longer 
periods have larger amplitudes. This ts the trme domiain analogy of the f-? relationship 
found in the spectra of bLODysno and OLOD,,,. The lowest period that correlation can 
be expected is around 20 days because the 23-point filter used to smooth the data 
completely removes periods of I5 days or less. Additionally, contamination from the 
unwanted tidal signal at a period of 13 days would make correlation at or near that pe- 
riod impossible. 
1. Seasonal Components 
Seasonal variations stand out, with the largest values of both OLOD,y..5 and 
OLOD,m Occurring in January and July. Munk and MacDonald (1960) observed these 
seasonal effects, and also noted that the amplitudes of the annual oscillations were 
greater in July than in January. This corresponds to the uneven contribution (by almost 
a factor of two in Munk and MacDonald’s computations) of the two hemispheres to the 
angular momentum budget during January and July. Thus, a typical day tn January ts 
almost one millisecond longer than a typical day in July. In January, strong westerlies 
from the winter (northern) hemisphere dominate the global angular momentum and 
cause the observed decrease in the length of day. In other words, the Earth spins slower 
because the atmosphere is spinning faster. In July, the upper level westerlies in the 
northern hemisphere are weakened by the breakdown of the midlatitude baroclinic zone. 
Meanwhile, the southern hemisphere, because of its smaller land mass, does not com- 
pensate with stronger westerlies. The net effect is for a minimum of atinospheric angular 
momentum in July and a corresponding maximum in Earth rotation rate, and the reverse 


in January; maximum angular momentum and minimum earth rotation rate. 


6 The linear correlation coefficient is defined as p,, = Cov,,/o,c,. Its analogy in frequency 
space is coherence, although the coherence statistic used in this ao is more Pe called co- 
herence squared (y 2); see Eq. (20). Direct comparisons of the two statistics, y ? and p, is of qual- 
itative value only because y? 1s frequency dependent whule p,, 1s not. 
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The seasonal components of 0LOD...,5 and 0LOD,,, (fig. 4 on page 23) appear 
to be two in-phase sinusoids of annual and semiannual frequency. The linear correlation. 
0.990, is significantly higher than the subseasonal linear correlation of 0.853 (refer to 
Table 1). The 0\LOD,..,5 has slightly greater amplitude, shown bv its greater variance 
(Table 1). Since a fundamental assumption (inade earlier) is that the OLOD,<.9 signal 
represents the total of all atmospheric (as well as other, unmodelled, sources of) angular 
Inomentum, it 1§ not surprising that the seasonal terms of the atmosphere are less than 
the tune signals. The neglected pressure term in (13), oceanic currents, and stratospheric 
winds are probably the missing inputs to balance the seasonal terms. Of particular note 
is the difference in amplitudes found in mid July of 1984 and 1986. The Quasi- 
Biennual-Oscillation (QBO) peaked at exactly those two times (R.A. Madden, personal 
communication, 1988). The peak was for easterly stratospheric winds; enough to 
qualitatively account for the significant amplitude difference of the troposphere and 
measured LOD. 

Figs. 5 and 6 are the decompositions of the seasonal signals depicted in Fig. 4 
on page 23 into annual and semiannual parts. Annual terms, Fig. 5 on page 24, appear 
smaller than those determined by harmonic analysis (Table 3) because the filter transfer 
function (required to separate the closely spaced annual and semiannual terms) atten- 
uated the annual signal 65%. Fig. 19 in Appendix C is the plot of the spectral response 
of the filter used to isolate these terms. Residuals of 6LOD,sx5 and dLOD,,,, show a 
small biannual bias which is probably associated with the QBO. 

Fig. 6 on page 24 contains the semiannual part of the seasonal signal. The 
spectral response of the filtering (Fig. 20, Appendix C) show minimal amplitude losses 
due to the filtering process. The amplitude from the plot, about 0.3 ms, agrees well with 
the harmonic analysis (Table 3) and with the findings of Rosen and Salstein (1985). 
These data show, as other data have in previous studies, that variations in the length of 
the day on seasonal timescales (annual and semiannual periods) are almost entirely due 
to meteorological influences. 

2. Subseasonal Components 

Subseasonal 1s defined to include periods from 20 to 150 davs. Fig. 7 on page 
25 shows the close match of the subseasonal signals in both phase and amplitude. The 
correlation coefficient of the two series is 0.853. The high peaks are associated with pe- 
riods around 50 days while the lower peaks correspond to shorter periods, as 1s expected 


with the f-? power spectrum discussed earlier. The 50-day peaks are perhaps associated 
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with the 40-50 day oscillations in the tropics first described by Madden and Julian 
(1971). The average amplitude of 0.25 ms ts slightly larger than that found by Eubanks 
et al. (1985) for the period 1977-1981, and agree well with the findings of Morgan et al. 
(1985) for the period 1981-1983. 


B. GLOBAL COMPARISONS IN THE FREQUENCY DOMAIN 
i. The Energy Density Spectra of OLOD,,,, and LOD, <x 

Plots of dLODy 59 and dLOD,,, energy density, coherence and phase are con- 
tained in Fig. 12 on page 28. The most striking feature of the energy density plot is the 
large peaks at annual and semiannual periods. This is consistent with the largest ampli- 
tude oscillations in the time domain (Fig. 2 on page 22) occurring seasonally. Also 
noteworthy are the sharp spikes at periods of 13 and seven days in the 0DLODy 5x9 Spec- 
trum. These are high frequency tidal oscillations which were not filtered out of the ori- 
ginal time data. 

The slope of both spectra is -2 (indicated bv the sloping double line on the en- 
ergy density plot). Eubanks et al. (1985), state that, while the cause of the f? dependence 
in the power spectrum 1s not clear, it shows up in a wide variety of geophysical phe- 
nomena. Near the Nyquist frequency of ‘Stay the OLOD Ln. spectrum falls off at a 
more rapid rate than expected by the f? relationship. This is an indication of the lower 
noise level of the 6LOD,.., data. Eubanks er al. (1985) experiments were conducted 
prior to recent advances in more accurate astronomic LOD measurements; thus their 
data probably have a somewhat noisier spectrum at higher frequencies (personal com- 
munication with J.O. Dickey, co-author with Eubanks ef al., 1985). 

At periods of 16 days and greater, the two spectra are in agreement. In addition 
to the broad peaks at annual and semiannual periods, there are broad peaks at 75 and 
50 days, and narrower peaks at 33 and 25 days. 

2. Coherence Spectrum 

The coherence spectrum in Fig. 12 shows strong coherence (or correlation) at 
seasonal as well as periods centered around 75, 50, 33, 25 and 20 days. As would be 
expected, these periods show corresponding power in the energy density plots. The 95% 
confidence line in the coherence plot is derived from a 95% significance level in rejecting 
a null hypothesis of no coherence between the two series. The coherence plot agrees 
well with that of Eubanks er al. (1985). Of note is the rapid drop off in coherence at 
periods around 125 days (also found by Eubanks et al. 1985). This may be related to a 
general lack of power in both the dLOD,,,, and 6LODysno spectra, however, the loss of 


power does not completely explain the steep drop in coherence. It is possible that 
whatever (weak) oscillations there are in the length of day at periods around 125 davs, 
they are not of meteorological origin. Significant coherence at periods as low as 20 days 
indicates that the atmosphere can generate short period, high frequency, changes in the 
LOD. Given accurate enough time and wind data, fluctuations of even one or two day 
time scales will eventually be resolved (Lambeck, 1980). 
3. Phase Spectrum 

The phase plot shows that, except for a small spike at 125 days, changes of 
6LOD,,,and dLODysy9 Occur Within five days of each other. Generally, phase is closest 
to zero Where coherence 1s highest, indicating changes in angular momentum are rapidly 
transmitted between Earth and atmosphere. Measurements of atmospheric angular 
momentum in this study cannot readily be decomposed into specific torques responsible 
for transfer of momentum. Swinbank (1985) found that shorter period fluctuations 
(from daily to subseasonal) result from mountain torques (pressure differences across 
north/south oriented mountain ranges) and seasonal variations result from friction 
torques (frictional drag at the surface boundary layer of the Earth and atmosphere). 

4. Harmonic Analysis 

Amplitude and phase information of the first two harmonic frequencies of 
OLODusxo and OLOD,,, were calculated by vear. The results (Table 3) compare favora- 
bly with previous studies (Morgan et al., 1985, and Rosen and Salstein, 1985), although 
the annual amplitudes found in this study show more variability than those found in the 
previously cited studies. Of note is the two-month variation in phases (Jan-Feb for an- 


nual terms and May-Jun for semiannual terms). 


C. LATITUDE COMPARISONS OF LOD 

Comparisons of 0dLOD,;x9 With tropical and extratropical partitions of 0LOD,,,,are 
shown in Figs. 13 through 16. As noted above, the reason for computing 06LOD,,,, over 
specific latitude bands (equatorial and extratropical regions, respectively) is to try to 
identify the particular geographical region that is contributing to the high coherence 
between OLOD,,,, and 0LODysyo. Energy density plots show significant drops in power 
at periods between one year and 30 days for the extratropical, 6LOD,,,75 and OLOD, a1; 
series. The two tropical series, (6LOD nis and OLOD, in.) Show power losses only at pe- 
riods longer than 100 days. Coherence for the extratropical series dropped below 0.5 for 
most periods between 30 and 100 days; for the tropical series, coherence at these periods 


dropped only slightly. Coherence falls off rapidly in both tropical and extratropical series 


1? 


around periods of 125 days. The extratropical series fall to near zero at these periods. 
This may be related to the local minimum of power at periods from 80 to 125 days in 
all atmospheric and USNO energy density plots. The phase plots show that at the an- 
nual period, changes in the atmosphere tend to lead changes in the length of day in 
extratropical latitudes (Fig. 14 and especially 16) while they tend to lag changes in the 
length of day in the equatorial regions. This simply reflects the fact that the annual cycle 
in the equatorial region 1s shifted or delayed somewhat in time relative to the middle 
latitudes. 

At subseasonal periods, it is very difficult to interpret much about the geographical 
origin of events with periods shorter than about 30 days. The generally high power and 
coherence near 20-30 days in Fig. 12 does not seem to be of only equatorial or only 
midlatitude origin. However, the evidence very clearly indicates that the tropics and 
equatorial regions are responsible for most of the fluctuations of dLOD in the 40-to 
{00-day range. Thus, the high power and coherence at 50-60 day and at 70-100 day 
periods in Fig. 12 also appear in Figs. 13 and 15 (using only tropical wind data), but they 
are completely absent in Figs. 14 and 16 where only extratropical wind data were used 
in computing dLOD,,,. The dividing line between the tropical and extratropical regions 
is not very clear because only two latitude bands have been considered. Also, seasonal 
signals are still apparent in the tropical regions, and, to a lesser degree, subseasonal 
(40-to 100-day periods) oscillations are seen outside the tropics. 
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IV. SUMMARY 


[t is clear that changes in atmospheric angular momentum can account for the vast 
majority of short period (one vear or less) changes tn the LOD . Time and frequency 
domain analvses show that over 90% of the variance in seasonal and subseasonal peri- 
ods of OLOD,,,, is explained bv 0LOD .y5. Meteorological activity 1s responsible for all 
changes in LOD within present measurement capabilities. While error analysis of the 
atmospheric and time series ts beyond the scope of this paper, Morgan er al. (1985) ac- 
counted for most of the variances between the series as being within measurement limits. 
In addition, Rosen and Salstein (1985) reduced the variance due to unknown sources to 
less than 0.04 (ms)? for annual and 0.03 (ms) for semiannual terms (measured tn units 
of amplitude). Phase differences between 6LOD,,,, and 0LODysyo are Very near zero, in- 
dicating that an almost instantaneous transfer of angular momentum occurs between 
Earth and the atmosphere. Additionally, when the pressure term is included, an rms er- 
ror of less than 0.25 (ms)? at seasonal periods results (Barnes ez al., 1983). Other studies 
(Eubanks er al., 1985 and 1986), conclude that measurement errors exceed the unex- 
plained variance of dLOD. 

The value of establishing the strong coherence and correlation of the independently 


measured signals of 6LOD,,,, and dLODysx9 include: 


e Astronomic measures of LOD can be used as an independent measure of the state 
of the zonal atmospheric winds. Analysis and forecasting models can use this as a 
check of data consistency. 


e Rapid changes of LOD, such as were observed in the El Nino event of 1983, can 
perhaps be used to infer future E] Nino events (Eubanks er al. 1986, and Rosen 
etal. 1984). Other weather phenomena with a strong zonal component, such as the 
40-50 day tropical oscillations, can be temporally located by real time filtering and 
analysis of 6LOD,,,, data. 


e Dickey ez al. (1986) reported that atmospheric angular momentum can be used to 
infer real trme OLOD. This is required for extremely precise Earth position fixes to 
navigate spacecraft. Real time calculations of atmospheric angular momentum and 
OLOD am are Currently computed and archived by NMC. 

The findings of this paper are in agreement with previous studies. The finer resol- 
ution of the USNO time data beginning in 1983 (the BIH time data also had similar 


improvements starting in 1983) shows that atmospherically and astronomically derived 
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OLOD are coherent down to periods of 20 days or less. [he NOGAPS wind analyses 
compare well with similar studies done with NMC or ECMWE. 

Comparing latitude bands of integrated 6LOD,,,, with bLOD¢esxo Confirms that there 
is spatial variability in both phase and amplitude of atmospheric angular momentum. 
Seasonal changes at annual periods of 6LOD appear to originate primarily in the 
extratropical regions. Subseasonal fluctuations (30-100 day periods) in dLOD,,,have 
significantly more amplitude (power) and coherence with dLOD,..9 when computed 
from tropical data than extratropical data. This strongly suggests that changes in the 
length of day at these time scales 1s produced by meteorological events in the tropics- 
e.g. the 30-60 day oscillations of Madden and Julian (1971, 1972). Further study, in- 
cluding more latitude band divisions of the globe is required to be more exact on the 
position of the dividing line between the tropics and extratropics. This study indicates 
that the latitude band responsible for most of the subseasonal fluctuations of the length 


of the day 1s approximately 20°N to 20°S. 


Table 1. LINEAR CORRELATION COEPMCIEM> 
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Table 2. VARIANCES 
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Table 3. PHASE AND AMPLITUDES OF LOD 


Annual | Semiannual 
of amplitude! | phase? | amplitude [phase 
May 12 
50a May7 
0OLOD ye 


Note 1. Amplitude 1s in ms. 

Note 2. Phase is the peak value of that term in days starting with Jan 1. 
For the semiannual phase, the date is the day of the first peak 
in the year. 

Note 3. From Rosen and Salstein (1985). Time data is from BIH. NMC 
weather data extends to | mb. 

Note 4 From Morgan er al. (1985). Time data is a composite of 
LLR, VLBI, and BIH estimates. NMC data extends to 100 mb. 
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Heavv line is USNO series. 


Seasonal terms of OLOD,,, and dLOD,ysyo: 
A 75-point weighted moving average filter was applied twice to each series; 


148 points were lost from the ends. 


Fig. 4 
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Heavy line 1s 0DLOD,..9>Bottom 


period of 2 years and 1s probably related 
Y 


to the Quasi-Biennual-Oscillation. Frequency transfer function is Fig. 19 


on page 48 of Appendix C. 
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Heavy line is 0LOD¢sno 


vo - OLOD im shifted down by 0.5 ms. 
0 on page 48 of Appendix C. 
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Semiannual terms of 0DLOD.<., and dLOD,,,' 


Bottom graph is residual of dLOD,, 


Frequency response of filter is Fig. 


Fig. 6. 
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Residual of series 


depicted in Fig. 4 on page 23 subtracted from Fig. 3 on page 23. Bottom 


Subseasonal Oscillations of OLOD,,, and dLOD:sxo0: 
line is the difference of the two series shifted down 0.5 ms. 


Fig. 7. 
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Heavy line is 0OLOD esxo 


seasonal signal: light line is OLOD,,,.;; derived from integrating atmospheric 


Comparison Of OLODysxo-subeees ANd OLOD amis: 
angular momentum within 15° of the equator. 


Fig. 8. 
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Comparison of DLODusxo-ccacona 2Nd OLOD yas? 


midiatitude components poleward of 1 
seasonal terms of OLOD 5x0. 
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NO sub- 


seasonal signal (from Fig. 7); light line is derived from integrating atmo- 


Heavy line is US 
spheric angular momentum within 20° of the equator. 


Comparison of 0OLOD esyosubsers 2d OLOD ymr0? 
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); light line is atmospheric signal derived from 


midlatitude components of atmospheric LOD above 20° VN and below 
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LOG 10 OF POWER SPECTRAL DENSITY 


COHERENCE 


OeZ5 


PHASE 


Fig. 
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1000 365 162. 160: 60 94050) * 205 is15 5 


Energy Density, Coherence, and Phase Spectra of OLOD,,, and 
SLODesxo: Seven point spectral smoothing used for all plots. Solid line 
in the power spectral density spectral plot is OLODysxo and dashed line 1s 
56LOD,,,. Dashed lines in phase plot are 5 day advances (top dashed line) 
and lags (bottom dashed line) of 6LOD¢sxo with respect to OLOD mn. 
Horizontal line in coherence plot represents 95% confidence that the two 
series are correlated. 
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13. Energy Density, Coherence, and Phase Spectra of dLOD,sx9 and 
OLOD hms: AS in Fig. 12 except the atmospheric signal is 0OLOD,,,, com- 
puted within 15° of the equator. 
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Fig. 15. 
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APPENDIX A. FORTRAN CODE FOR DECODING NEDN TAPES AND 


COMPUTING LOD 


//BENEDICT JOB (4507,0818),'JAN 87',CLASS=G 
//*MAIN ORG=NPGVM1. 4507P, LINES=20 
// EXEC FORTVCLG, PARM. FORT='OPT(2),LVL(66)',REGION. GO=2048K 


PeemORT.SYSIN DD * 
Cede ts dove Toe ve He ve ve Te ok ae 96 Fe He Ie Me He Vere Ieee TM POR TAN T 20 ae 3k eae ve 26 Fe 90 ae ae Te se eae He ak ae ae He He Ie Ie ae oe Ve He He He 


C 


Co G2 C2 Cl GaGa C2 C2 C2 C2 CD 


ORIGINAL AUTHOR JAMES BOYLE, NAVAL POSTGRADUATE SCHOOL 1986 


ADDITIONS MADE IN 1987 BY W. L. BENEDICT 


THIS 


AVERAGED WINDS OUT. 
END OF THE PROGRAM FOR THE DATES. 


IS THE MAIN PROGRAM TO READ THE NEDN TAPES AND GET ZONALLY 
BE SURE TO CHECK THE INPUT PARAMETERS AT THE 
THE OUTPUT FILE WILL GO TO YOUR 


READER. 
FeV TTT TTT ETAT RTE EEEEERERE TMPORTAN T reve tere ce tere rete cede Fede seve cede eT Tee Te Teh RRC 


PGM TO READ DATA FROM FNOC NOGAPS TAPE GLOBAL GRID NEDN FORMAT 


>>>>>>>>>>>>>>>>>>> LONG RECORD FORMAT <<<<<<<<<<<<<eceeeeeeccc 


REAL*4 DATAC 10512) ,DC73,144) ,C( 144,73) 
REAL*4 ZONAVG( 11,84) ,COSSQ( 73) ,AVERGE , FACTRV(8) , VERT( 73), 


& 


LODATM(400), OMEGA, ISHELL,AVG,G,PI,LODSID,R 


INTEGER IBUF(5284) 

INTEGER IDATA( 10512) 

INTEGER*2 JDATA( 21024) 

INTEGER IN1( 24) ,KFTC( 100) 

INTEGER BCD( 24) ,FTC( 100) 

IVIEGER LOPELDC 3,12) ,KOUNIC6,12,6) 
INTEGER DELTAU, COUNT ,QQ,Q 

INTEGER CALNDR( 740) 

REAL FCTR( 12) ,DBASE( 12) 

INTEGER ITAUOK(4) 


EQUIVALENCE (IDATA(1),JDATA(1)) 
EQUIVALENCE (DATA(1),D(1,1)) 


NAMELIST/INPUT/ NSKPRC,IX1,1X2,IDF1,IDF2 


DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 


FCTR/ 12%0.01/ 
DBASE/12*0. 0/ 
ITAUOK/00,12,24,36/ 
NUMTAU/01/ 
DELTAU/12/ 

FACTRV /2%7500, 10000, 3*5000, 2*2500/ 
OMEGA /7. 29E-5/ 
ISHELL /7. 04E37/ 

PI /3.1415926/ 

G /9.81/ 

LODSID /86000. 0/ 
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C 


iy 


1276 


DATA R /6370000. 0/ 


READ( 05, INPUT) 
WRITE( 06, INPUT) 


READ(0S,1275) (CCIDFELD(1,J) =o ote DE ph?) 
FORMAT(11(3A1,1X)) 

WRITE( 06,1276) (CIDFELD(1,J) , 1=1,3).0=I1DEP Loe) 
FORMAT(1X,'FIELDS SELECTED ',11(3A1,1X)) 


COUNT=0 
QQ=0 
NA=21135 
NB=0 
NCNT=0 
NF LE=0 


C SET UP ARRAY FOR CALENDAR 


CALL KLNDRCCALNDR) 


C SET UP COSINE ARRAY FOR INTEGRATION 


C 


7 


BS LG TD I Se I cA 


3G) 


2300 


10 


DO 2300 I=1, 37 
COSSQ(I) = (SIN(2. 5%*(I-1)*(2*PI/360. 0) )**2) 
COSSQ(74-1) = COSSQ(I) 

CONTINUE 


DO 10 I=1,6 
DO 10 J=1,8 
DO 10 K=1,6 
KOUNT(I,J,K)=0 


IDDHH1=IX1 - ((1IX1/10000)*10000 ) 
TIDD1=IDDHH1/100 
THH1=IDDHH1 - IDD1*100 


CALL BUFFER(10,1,21135, IBUF) 


SKIP TO DESIRED DATA E.G. DEC 74 NSKPRC = 806 


CALL SKPREC(10,NSKPRC,&99 ,&999) 


IDFELD IS THE NEDN FIELD IDENTIFIER E.G. AOQl = SLP 
FQ0=500 Z, T21 = 250MB V, 


IX] 


& IX2 ARE THE DTG OF FIRST AND LAST PERIODS (FORM YYMMDDHH IN 1[8) 


Week: TOP OF DATA INPUT LOOP vericidedetietedcikk 


55 


CONTINUE 


LONG RECORD FORMAT JUST USE ONE READ 


CALL GETREC(10,&99 ,&9999) 
NCNT=NCNT + 1 


CALL GBYTES( IBUF(1),IN1(1),00,06,00, 24) 
CALL XBCD(24,1IN1(1),BCD(1)) 

CALL GBYTES( IBUF(6) , IFTC, 28,4,0,0) 
NFTC=IFTC*#4 | 
CALL GBYTES( IBUF(16),KFTC(1),24,6,0,NFTC) 
CALL XBCD(NFTC,KFTC(1),FTC(1)) 
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3G" C2 


C)Ga Co Ci-G) 


CALL XINT(BCD,5,4,ITAU) 


seaicskescveseocvevevesdeore FORECAST HOUR ( TAU ) CHECK Yededevesvasde vesevevedevesacdeskc deve de ve ve sc oF ve 


DO 610 K=1,NUMTAU 
NNTAU=K 
IFC ITAU. EQ. ITAUOK(K)) GOTO 620 


610 CONTINUE 


620 


GOTO 55 
CONTINUE 


weve’ =CHECK YYMMDDHH TO GET THE DATA WE NEED = *%%o%% 


CALL XINT(BCD,9,8, IT YMDH) 

PEG. LIStAl) GO TO 55 

eel: (nee x2) GO TQO./7 

IDDHHX=IYMDH - (CIYMDH/10000)*10000) 

IDDX=I DDHHX/ 100 

THHX=IDDHHX - IDDX*100 

NNDAY=(IDDX-IDD1) + CIABS( IHHX-IHH1)/12) + 1 


‘eloreiedekeye DATA TYPE CHECK Yevedevedededededesevevevevdevesedc doves vex 


644 
640 
645 


desese 


Yeveve 


CDC 


30 


DO 640 K=IDF1,IDF2 

NNFLD=K 

DO 644 KK=1,3 

IF(BCD(KK) .NE. IDFELD(KK,K)) GOTO 640 
GOTO 645 

CONTINUE 

GOTO 55 

CONTINUE 


DECODE SCALING INFORMATION, GRID SIZE 

CALL GBYTES( IBUF(6),IUNIT,8,5,0,0) 

CALL GBYTES( IBUF(6),ISCLE,13,6,0,0) 

CALL GBYTES( IBUF(6) ,KEYBIT, 26,2,0,0) 

CALL GBYTES( IBUF(7),M,24,12,0,0) 

CALL GBYTES(IBUF(8) ,N,4,12,0,0) 

NR=M*N 

DECODE DATA ITSELF AFTER FINDING WHERE IT STARTS. 
NBITS = (21 + IFTC)*24 

IWDSRT = (NBITS/32) + 1 

ISKIP = NBITS - ((NBITS/32)*32) 

NBTSDP=16 

CALL GBYTES( IBUF( IWDSRT) , IDATA, ISKIP,NBTSDP,0,NR) 
CALL GBYTES( IBUF(27) , IDATA,8,NBTSDP,0,NR) 


HAS ONE'S COMPLEMENT, IBM HAS TWO'S COMPLEMENT FOR NEGATIVE 'S. 
DO 30 I = 1,NR 
IF (IDATACI) .LT. 0) IDATACI) = IDATACI) + 1 


SET UP TO UNPACK 16 BIT DATA POINTS. THUS JDATA WHICH IS *2 HAS ONE 
DATA POINT PER WORD. DATA IS IN THE RIGHTMOST 16 BITS OF IDATA WHICH 
IS 4, 


L=0 
NR2=NR*2 
IF (ISCLE .LT. 32) GO TO 40 
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ISL= 2**(1Seumeg2) 
DO 17 1 = Seno 
L=L+1 
IDUM=JDATA(I) 
17. DATA(L) = FLOAT( IDUM) /FLOAT( ISL) 
GO TO 45 
40 ISL=2%*( 32-ISCLE) 
DO 18 I=2,NR2,2 
L=L+1 
IDUM=JDATA(I) 
18 DATA(L)=FLOAT( IDUM)*FLOAT( ISL) 
45 CONTINUE 


C 
C DBASE = 5574.0 FCTR = 0.01 FOR SOOMB FOO 
C DBASE = 762.0 FCTR = 0.01 FOR 925MB ROO 
C DBASE = 0000.0 FCTR = 1.00 FOR SLP AO1 
C DBASE = 0000.0 FCTR =+0.01 FOR V AND 0.01 FOR JU. 
DO 88 KZ=1,NR 
88 DATACKZ)= DATACKZ)*FCTRCNNFLD) + DBASEC(NNFLD) 
C 
C ORIENT I,J INDICES SUCH THAT I=E-W J=N-S 
C REVERSE SENSE OF J SOQ THAT IT INCREASES NWD. 


OZomi—lc 
II=73 - (I-1) 
DO 25 J=1,144 
256 Cl is = UGiee) 


IHR=IYMDH-( IYMDH/100)*100 
IDAY=( ( LYMDH-( IYMDH/ 10000)*10000) -IHR)/100 
C  werteredededededeteiedcdick ZONAL AVERAGING LOOP xe vetevededeacde te dete edeve kde ves seve ve tere eH He HH He He 
C 
2999 DO 2a) = lies 
AVG = 0.0 
DO 2116 I= 1, 144 
AVG = C(I,J) + AVG 
2116 CONTINUE 
ZONAVG(NNFLD,J) = AVG/144. 
Z2liS CON ANGE 
CG tetedededededededesededee de ve ae te ve te Fe te se ve ve He se ve He He He He Fe He He He Fe He He He Fe He He He Fe He He Fe Fe He Ke Fe Fe Ke Fe Fe Fe He Fe Fe He Ve He He Me Fee He 
IF (NNFLD .EQ. 11) GOTO 3000 
GOTO 55 
3000 QQ = QQ +1 
IF (CIYMDH .NE. CALNDR(QQ)) GOTO 3500 
GOTO 3600 
3500 WRITE(6,3501) QQ, CALNDR( QQ) , LODATM( QQ) 
3501 FORMAT(1X,13,1X,1I8,1X,F9.6,2X,'DATA MISSING FROM FNOC TAPES' ) 
GOTO 3000 
C wevededededodedededededotktedekickeRe VERTICALLY INTEGRATE %% vtec re vetedececeve te eve seve dete secede deserve 
3600 DO 2130 I=1, 73 
VERT(I) = (ZONAVGC1,1I)+ZONAVG( 2,1) )*FACTRV(C1) 
( ZONAVG( 2, 1)+ZONAVG( 3,1) )*FACTRV(2) 
(ZONAVG( 3, 1)+ZONAVG(4, 1) )*FACTRV(3) 
(ZONAVG(4 ,I)+ZONAVG(5 , I) )*FACTRV(4) 
(ZONAVG(5 , I1)+ZONAVG( 6,1) )*FACTRV(5) 
( ZONAVG( 6 ,1I)+ZONAVG(7 , I) )*FACTRV(6) 
(ZONAVG( 7, I} +ZONAVG( 8,1) )*FACTRV( 7) 


MOO OEE 
tetttt 


36 


& + (ZONAVG(8,1I)+ZONAVG(9,1))*FACTRV( 8) 
& + (ZONAVG(7,1)+ZONAVG(8,1I))*FACTRV(9) 
2130 CONTINUE 


C APPLY COS LAT TERMS 
DO 2140 I=1, 73 
2140 VERT(I) = VERT(1I)*COSSQ(I) 


c HORIZONTALLY INTEGRATE 
AVERGE=0. 0 
DO 2145 I=2,72 
AVERGE = VERT(I) + AVERGE 
2145 CONTINUE 
AVERGE = AVERGE* PI/72.0 
C 
C CONSTANTS OUTSIDE THE INTEGRAL ARE NOW BROUGHT IN: 
AVERGE = AVERGE*(2*PI*(R**3))/G 
C CALCULATE THE CHANGE IN LOD 
LODATM(QQ) = 1000. *(AVERGE*LODSID/( OMEGA*ISHELL) ) 
C 
WRITE(6,3001) QQ, IYMDH, LODATM(QQ) 
3001 FORMAT(1X,1I3,1X,1I8,1X,F9. 6) 
GOTO 55 


99 CONTINUE 


PewRUTE( 6,151) LYMDH 
CC vkevederedededeverevevevereoedcderedererek QUTPUT LOOP wktedededededcsedkeincikkkekikkkkeKee 


WRITE (6,1020) 
C WRITE (6,1010) (LODATM(Q) ,Q=1,QQ) 
1010 FORMAT(7F10. 6) 
1020 FORMAT(1X,'THESE ARE THE FINAL LOD VALUES IN SECONDS' ) 
151 FORMAT(1X,'PGM ENDED OK IYMDH = ',1I10) 


DO 835 I=1,NNDAY 

WRITE( 06,1235) I,1X1 
1235 FORMAT(1X,//,1X, NDAY = ‘,13,2X, ‘1X1 = ',I8) 
835 CONTINUE 


C2 C-o2 


CALL FILEND( IDF1, IDF2) 
STOP 
G Ks 
CG vededetstcdedsvedeedede dete PARITY ERROR TRAP %reverededededededededededededevedevededeck tev 
9999 WRITE(06,120) 
120 FORMAT(1X,' PARITY ERROR GETREC ‘) 
STOP 
G Sevevevevevedecedseteseseteseveve SKPREC FAILURE yevededeteteseteteseddestcksekedtedete tk vesede 
999 WRITE(06,122) 
122 FORMAT(1X,' SKPREC FAILURE ') 
STOP 
998 CONTINUE 
WRITE( 06,1298) NFLE 
1298 FORMAT(1X,' PGM ENDED. NUMBER OF FILES PROCESSED = ‘,I4) 
C CALL FILEND( IDF1,IDF2) 
STOP 


3H 


END 
C &&&&SESE&EEEEEEEES&& SUBROUTINE XINT S&&&&6&6&&&6665565& 
C 
SUBROUTINE XINT(BCD,ISRT,NUM, IOUT) 
C THIS ROUTINE TAKES NUM EBDIC CHARACTERS FROM BCD 
C STARTING AT ISRT AND CONVERTS IT INTO A INTEGER NUMBER PUT INTO IOUT. 
INTEGER BCD( 24) ,DIGITS(11) 
INTEGER ey oie at eaee 
DATA DIGITS/ 0','1 ,'2','3', 4', 55 6) eeuenes: oe nee 
DATA INTGRS/O, 7 235450057 SOR 
DATA IBLANK/' yy, 


IEND=ISRT + NUM - 1 
DO 100 I=ISRT,IEND 
II =I - ISRT+1 
DO 200 J=1,11 
IF(BCD(I).NE. DIGITS(J))GO TO 200 
INTSCII)=INTGRS( J) 
GO TO 100 
200 CONTINUE 
Ec WRITE(6, 1000)II, ISRT,NUM, IEND,( BCD(K) ,K=ISRT, IEND) 
C1000 FORMAT(2H ,'ERROR,CHECK YOUR SUBROUTINE XINT' ,4(2X,14),2X,24A1) 
C STOP 
INTS(II) = 
100 CONTINUE 
LOUT=00 
IFCTR=1 
DO 300 I=1,NUM 
II = NUM -I+1 
IOUT = IOUT + INTS(II)*IFCTR 
IFCTR = IFCTR*10 
300 CONTINUE 
RETURN 
END 
C &&&&E&ES&EEESEEEEEEEE& SUBROUTINE XBCD &&&&&&&&&&&&&&& SEE EEEEEEEEEEEE 
E 
SUBROUTINE XBCD(N, INPUT, BCDOUT) 


SUBROUTINE XBCD CONVERTS N INPUT CHARACTERS FROM CDC EXTERNAL BCD 
TO IBM EBCDIC. THE INPUT CHARACTERS ARE STORED IN ARRAY INPUT 
AND THE OUTPUT (CONVERTED) CHARACTERS ARE STORED IN ARRAY BCDOUT. 


GC) C2 C2 © 


INTEGER HEX(47)/231, 232,233,234, 235,236, 237,238, 239,221,222,223, 
* 224,225 ,226,227,228,229,012,213,414,415,216,21/7,218,219 , 200 
* Z02,203,204,Z05,Z06,Z07,2Z08,2Z09,230,220,210,211,23B, 228,220 90 
* Z1C,Z0B,Z0C/, INPUT(N) 

REAL CHAR(47)/'A','B','C','D','E','F' 'G! | H 


Edie: ae ke Ee 
ye 'y! N Oe le ig! Re tot ir ty! iy! Wi, ryt 7 tot ' ? hy! 
x Pe 5 ee a Og tt gt ngs ati ont et ee roth tat ryt 
3 9 ’ 3 3 3 3 3 
og , #','@' fk, BCDOUT(N) 
ms "100 I=1,N 


DO 200 J=1,47 
IF( INPUT(I). EQ. HEX(J)) GO TO 150 
200 CONTINUE 
J=39 
150 BCDOUT( I1)=CHAR(J) 
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100 CONTINUE 
RETURN 
END 
ic Wired KKK KEK SUBROUTINE KLNDR wrirdededersdedcceve de ded teas deve de teak vere veut Hee ve ve 1 
SUBROUTINE KLNDRCCALNDR) 
INTEGER CALNDR( 740) ,0 
CG JAN (31 DAYS) 
CALNDR( 1) = 87010100 
DO 5001 I[=2,62,2 
CALNDRCI)=CALNDRCI-1) + 12 
CALNDRC I+1)=CALNDRCI-1) + 100 
5001 CONTINUE 
meerehd (28 DAYS) 
CALNDR(63) = 87020100 
DO 5002 [=64,118,2 
CALNDRCI)=CALNDR(I-1) + 12 
CALNDRC I+1)=CALNDR(CI-1) + 100 
5002 CONTINUE 
C MAR (31 DAYS) 
CALNDR(119) = 87030100 
DO 5003 I=120,180,2 
CALNDRCI)=CALNDR(I-1) + 12 
CALNDR( I+1)=CALNDR(I-1) + 100 
5003 CONTINUE 
C APR (30 DAYS) 
CALNDR( 181) = 87040100 
DO 5004 [=182,240,2 
CALNDR(C I)=CALNDR(I-1) + 12 
CALNDR(C I+1)=CALNDRCI-1) + 100 
5004 CONTINUE 
CGC MAY (31 DAYS) 
CALNDR( 241) = 87050100 
DO 5005 [=242,302,2 
CALNDR(I)=CALNDR(I-1) + 12 
CALNDR(C I+1)=CALNDR(I-1) + 100 
5005 CONTINUE 
Gee oUNe(30 DAYS) 
CALNDR(303) = 87060100 
DO 5006 I=304,362,2 
CALNDRCI)=CALNDR(I-1) + 12 
CALNDR( I+1)=CALNDR(I-1) + 100 
5006 CONTINUE 
C JUL (31 DAYS) 
CALNDR( 363) = 87070100 
DO 5007 [=364,424,2 
CALNDRCI)=CALNDR(I-1) + 12 
- CALNDR(I+1)=CALNDR(I-1) + 100 
5007 CONTINUE 
C AUG (31 DAYS) 
CALNDR( 425) = 87080100 
DO 5008 I[=426,486,2 
CALNDRC I)=CALNDR(I-1) + 12 
CALNDR( I+1)=CALNDR(I-1) + 100 
5008 CONTINUE 
C SEP (30 DAYS) 
CALNDR( 487) = 87090100 
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DO 5009 I=488,546,2 
CALNDR(CI)=CALNDR(I-1) + 12 
CALNDR(I+1)=CALNDR(I-1) + 100 

5009 CONTINUE 
C OCT (31 DAYS) 

CALNDR( 547) = 87100100 

DO 5010 I[=548 ,608,2 
CALNDRCI)=CALNDR(I-1) + 12 
CALNDR( I+1)=CALNDR(I-1) + 100 

5010 CONTINUE 
C NOV (30 DAYS) 

CALNDR( 609) = 87110100 

DO 5011 I=610,668,2 
CALNDRCI)=CALNDR(I-1) + 12 
CALNDRC I+1)=CALNDR(I-1) + 100 

5011 CONTINUE 
C DECRG Das 

CALNDR( 669) = 87120100 

DO 5012 I=670,730,2 
CALNDRCI)=CALNDRCI-1) + 12 
CALNDR( I+1)=CALNDR(I-1) + 100 

5012 CONTINUE 
RETURN 
END 
C reat reredererevedesedereakevekiedereeekeicseHe SUBROUTINE FILEND *iriitddedcdictetedteaaiekicded 
SUBROUTINE FILEND(C IDSRT, IDEND) 
KK1=C( IDSRT + 1)*10 
KK2=( IDEND + 1)*10 
DO 778 KK=KK1,KK2,10 
DO 778 K=2,02 
IFILE = KK + (K = 1) 
778 ENDFILE IFILE 
ENDFILE 21 
ENDFILE 22 
ENDFILE 24 
ENDFILE 26 
RETURN 
END 


GG Gas) Gy). C2 C9 0) 


~ 


/LKED. SYSLIB DD 

/ DD 

/ DD 

Q DD 

/ DD 

i DD UNIT=3350 , VOL=SER=MVS003 , DISP=SHR , DSN=C1012. TAPEUT, 

/ LABEL=(,, , IN) 

//* NOTE ONLY ODD NUMBERED TAPES HAVE THE 00Z DATA AND 120HR FCSTS. 
//GO.FT10F001 DD UNIT=3400-5 , VOL=SER=W16733, 

Ta DISP=(OLD, PASS) , LABEL=(1,NL,,IN), 

// DCB=( RECFM=FB , BLKSIZE=21135 , DEN=4) 

/ (GO. SYSiN pie 

& INPUT NSKPRC=00, IX1=87070100, IX2=87123112, IDF1=01, IDF2=9&END 
C20 D20 E20 F20 G20 H20 I20 J20 K20 


//* 


~~, “A, TA, TA, TM, TM, 


4() 


APPENDIX B. FORTRAN CODE FOR SPECTRAL ANALYSIS 


Crete reverted Rie te Hs Te eK Ts Te Le Fede KC ETT T ETRE T ETE PERE LO TELE TE TET EET ETE TE VED EE TENTED TEEN TEBE TER TRIED EIR EER ICE 


C PieirGy VENollY 2s CRGss SEBCTRUM] PHASE, AND COHERENCE 

C PweCP FILES ARE UNITS 33, 34, and 35 

C OUTPUT FILES ARE UNITS 41 THROUGH 45 

C 

Cisdevededevedetessdedeceteds Lede Te LR Te PRT LTE TEE TCT FTE TS TET TC TE TENET TEC TED TE DET ED TE VE TE FEET CTE POTS FOREVER EET ER 
C 


PROGRAM CRSPEC 
PARAMETER( N=2048 , M=11) 


C 
C DECLARE ARRAYS. 
REAL*4 DATAX(N), GX(N), GXSUM(N/2 + 1),ATM1( 1622) ,ATM2( 1622) 
REAL*4 DATAY(N), GY(N), GYSUM(N/2 + 1),LOD( 1622), LIMPH1(N/2+1) 
REAL*4 DAY, XSUM(N/2+1), YSUM(N/2+1) 
REAL*4 CXY(N), QXY(N), PHASE(N), COHER(N), FREQ(N),LIMPH2(N/2+1) 
INTEGER*4 IWK(N+1) ,FILENO 
COMPLEX AX(N), AY(N), GXY(N), GXYSUM(N/2 + 1), XYSUM(N/2+1) 
Ee . 
CHARACTER*60 TITLE1,TITLE2, TITLE3,TITLE4, TITLES , TITLE6 , TITLE? 
C DEFINE CONSTANTS. 
C M AND N ARE DEFINED IN PARAMETER SECTION; N IS NUMBER OF SAMPLES 
C M IS THE POWER OF TWO TO GET N 
XN = FLOAT(N) 
NF = N/2 + 1 
C (NF IS INDEX VALUE FOR NYQUIST FREQ. ) 
pie 1 
C (DT IS TIME BETWEEN SAMPLED VALUES) 
T = FLOAT(N) 
c (T IS TOTAL TIME FOR SAMPLE; HERE IT IS EQUAL TO THE SAMPLE) 
DF = 1./T 
C (DF IS DELTA F, FUNDAMENTAL FREQ. ) 
NS = 1 
C (NS IS NO. OF SAMPLES TO ENSEMBLE AVERAGE) 
PI = 3.141593 
WCORR = . 875 
C (WINDOW CORRECTION) 
C 
C INITIALIZE SUMMING ARRAYS. RETAIN INFORMATION UP TO THE 
C NYQUIST FREQ. ONLY. 
DO 100 I = 1, NF 
GXSUM(I) = 0. 
GYSUM(I) = 0. 
Gxacttian) = (0. ,0.) 
100 CONTINUE 
C 


Crevevedededsdetededededede dese dete dete te tedede MAIN LOOP OF PROGRAM Sevevetedestetededtetek destetetetetetesetet teste ttek 


C COMPUTE ENERGY DENSITY SPECTRA AND ENERGY DENSITY CROSS SPECTRUM 
C FOR EACH SAMPLE AND ENSEMBLE AVERAGE. 
C 

FILENO=45 


41 


DO 500 K = 1, NS 


E 
C READ DATA. 
C THIS IS FOR 15N/S DATA 
DO 144 I= 1, 1622 
C READ (34,209) ATM1(I) 
C BELOW IS ONLY FOR 20N TO 20S DATA. 
C READ (33,209) ATM1(T) 
C ATM1(1)=16. /72. *ATM1(T) 
READ (35,200) LOD(I),ATM1(I) 
C ATM1(1)=ATM2(1)-ATM1(1I) 
144 CONTINUE 
C 


200 FORMAT (2(1X,F9.6)) 
202 FORMAT (24X,F8. 5) 
209 FORMAT(10X,F9. 6) 


C TAPER BOTH ENDS OF THE DATA 


CALL TAPER(1622,LOD) 

CALL TAPER( 1622 ,ATM1) 
C ADD ZEROS 

DO 152 I=1,N 

IF (C(I . LE. 213) Sener. Gl ssa Hen 
DATAX( 1I)=0. 
DATAY(I)=0. 

ELSE 
DATAX( I )=LOD( 1-213) 
DATAY( I )=ATM1( 1-213) 

ENDIF 

152 CONTINUE 


WRITE (37,101) (DATAY(I),I=1,N) 
101 FORMAT( 8F8. 4) 
TRANSFER THE DATA TO COMPLEX ARRAYS AX AND AY FOR FFT2C INPUT. 
NOTE THAT IN USING THE 2ND FORMULA IN FFT2C DOCUMENTATION, YOU 
MUST INITIALLY TAKE THE COMPLEX CONJUGATE OF THE DATA. SINCE 
THE IMAGINARY PART IS ZERO, THIS STEP IS NOT NECESSARY. 
DO 300 I =1, N 
AX(I) = CMPLX (DATAX(I), 0.) 
AY(I) = CMPLX (DATAY(I), 0. ) 
300 CONTINUE 


to (2 C212 


C 
C DO FFT. 
CALL FFT2C (AX, M, IWK) 
CALL FFT2C (AY, M, IWK) 
C 
C MAKE ADJUSTMENT TO GET 2ND FORMULA. 


DO .325015= TN 

AX(I) = CONJG (AX(I)) / XN 

AY(I) = CONJG (AY(I)) / XN 
325 CONTINUE 


COMPUTE THE ENERGY DENSITY SPECTRA GX(F) AND GY(F) AND THE 
CROSS SPECTRUM GXY(F) FROM FFT2C OUTPUT AX AND AY. 


2 | 
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Dome50 I = 1,.N 
Gx 1) CABS( AX( I) **2)*( 2. /DF) 
Crcl ) CABSCAY( I) **2)*( 2. /DF) 
GXYC(I) = AYCI)*CONJGCAX(CI))*(2. /DF) 


350 CONTINUE 


S36) GC) 


400 
C 
500 


SINCE YOU MUST ENSEMBLE AVERAGE THE RESULTS OVER THE NO. OF SAMPLES, 
ADD GX(F), GYCF), GXY(F) INTO SUMMING ARRAYS. NOTE THAT EACH 
FREQUENCY IS SUMMED INDEPENDENTLY. 


DO 400 I = 1, NF 

GXSUM(CI) = GXSUMCI) + GX(I) 
GYSUN(CI) = GYSUMCI) + GYCTI) 
GXYSUMCI) = GXYSUMCI) + GXY(I) 
CONTINUE 

CONTINUE 


Crereretedodederedevere vedere se reve ae He ve He END OF MAIN LOOP Fore Tore te Te ve ve se te aie ve He vere Fe ve ve ve Fe te te Te v'e He te te verve ve Feve 


Cac) C2) Cre 
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600 


AFTER SUMMING IN ALL OF THE ESTIMATES, DIVIDE BY THE NO. OF 
SAMPLES NS AND BY THE WINDOW CORRECTION FACTOR WCORR. 
COU ine CO=sEECiRUM CXY, THE QUAD-SPECTRUM QXY, THE 
PHASE SPECTRUM, AND THE COHERENCE SPECTRUM. CONVERT THE 
PHASE SPECTRUM TO DEGREES FOR PLOTTING. 


XNS = FLOAT(NS) 


MULTIPLY BY F%**2 


DO 600 I = 2, NF 


GXSUM(I) = (GXSUM(I) / (CXNS * WCORR))*((FLOAT( I-1)*DF)**2) 
GYSUM(I) = (GYSUMCI) / CXNS * WCORR))*( (FLOATCI-1)*DF)**2) 
GXYSUM(I) = (GXYSUM(I) / (XNS * WCORR))*((FLOAT( I-1)*DF)**2) 


CONTINUE 


C APPLY A 5 POINT SMOOTHING CURVE) 


O10 


DO 610 I[=3,NF-2 

XSUM( I)=(GXSUM( I-1)+GXSUM( I )+GXSUM( I+1)+GXSUM( I-2)+GXSUM( I+2))/5. 

YSUM( I)=(GYSUM( 1I-1)+GYSUM( I )+GYSUM( I+1)+GYSUM( 1-2)+GYSUMC I+2))/5. 

XYSUM( I)=(GXYSUM( I-1)+GXYSUM( I)+GXYSUM( I +1)+GXYSUM(I-2) 
TeMoonC L Aa) yaar 

CONTINUE 

XSUM( 2 )=(GXSUM( 2 )+GXSUM(3))/2. 

YSUM( 2)=(GYSUM( 2)+GYSUM( 3) ) /2. 

XYSUM( 2 )=(GXYSUM( 2 )+GXYSUM(3) )/2. 


XSUM( 3)=( GXSUM( 2 )+GXSUM( 3 )+GXSUM(4))/3. 
YSUM(3)=(GYSUM( 2 )+GYSUM(3)+GYSUM(4))/3. 
XYSUM(3)=(GXYSUM( 2 )+GXYSUM(3)+GXYSUM(4) )/3. 


C RETURN VALUES TO ARRAYS 


676 
C 


DO 676 I=1,NF 
GXSUMCI)=XSUM(T) 
GYSUM( I)=YSUMCT) 
GXYSUM( I )=XYSUM(T) 
CONTINUE 


43 


C CREATE FREQUENCY ARRAY AS MULTIPLES OF THE FUNDAMENTAL FREQ. DF. 
C NOTE THAT THE FIRST BIN IS THE MEAN AND THE FREQ. IS ZERO. 
DO 700 1 = ies 
FREQ(I) = FLOAT(I-1) * DF 
700 CONTINUE 
DO 615 I=2,NF 
CXY(1) = REAL(GXYSUM(I)) 
QXY(I) = AIMAG(GXYSUM(I)) 
IF (GXYC1) See. 0.) e THEN 
PHASE(I) = O. 
GOTO 679 
ENDIF 
PHASE(I) = ATAN2(QXY(I),CXY(1))*( 180. /PI) 
IF (PHASE(I) .GE. 90.) PHASE(1I)=90. 
IF (PHASE(I) . LE. -90.) PHASE(1)=-90. 
679 IF ((GXSUM(I)*GYSUM(I)) .EQ. 0.) THEN 
COHER( I )=0. 
GOTO 615 
ENDIF 
COHER(I) = ((CXY(1)**2)+( QXY(1)**2))/(GXSUM( 1)*GYSUM(I)) 
615 CONTINUE 
C DIVIDE BY F**2 
BO 6012 NE 
GXSUM( 1)=GXSUM(1)/( ( FLOAT( I-1)*DF)**2) 
GYSUM( I )=GYSUM(1I)/( ( FLOAT( I-1)*DF)**2) 
GXYSUM( I )=GXYSUM( 1)/( ( FLOAT( I-1)*DF)**2) 
611 CONTINUE 
C 
C SET LIMITS FOR PHASE PLOT 
DAY=5. 0 
DO 833 I=2,NF 
LIMPH1( I )=FREQ( I )*DAY*360 
IF (LIMPH1(1I).GE. 90.) LIMPH1(1)=90. 
LIMPH2(I)= -LIMPH1(I1) 
833 CONTINUE 
C PRINT OUT FREQ. , GXSUM(F), GYSUM(F), CXY(F), QXY(F), PHASE(F), 
C  COHER(F). 
1005 FORMAT(14X,'X' ,9X,'Y' ,9X,'CO-' ,6X, ‘QUAD’ ) 
1010 FORMAT(1X,' FREQUENCY’ ,2X,'SPECTRUM' ,3X,'SPECTRUM' ,2X,'SPECTRUM'’ , 
&2X,'SPECTRUM' ,4X,'PHASE' , 2X, ' COHERENCE’ ) 
1015 FORMAT(1X,72('-')) 
Ee WRITE( 36,1005) 
C WRITE( 36,1010) 
cS WRITE( 36,1015) 
DOesOC ie —ecene 
WRITE(FILENO,810)FREQ(I),GXSUM( 1) ,GYSUM( I) ,PHASE(TI), 
& COHER( 1) ,LIMPH1(1),LIMPH2(I) 
800 CONTINUE 
DO 802 I = 5, NF-5, 5 
DO 844 K=1,5 
WRITE( FILENO,810)FREQ(K+I1) ,GXSUM( I+3) ,GYSUM( I+3) , PHASE( I+3), 
& COHER( I+3) , LIMPH1(K+I), LIMPH2( K+I1) 
844 CONTINUE 
802 CONTINUE 
810 FORMAT(7F10.5) 


ren 


C 
C 


COMPUTE THE VARIANCE BY INTEGRATING GXSUM(F) AND GYSUMCF) 
NUMERICALLY BY SUMMING GXSUM(F)*DF AND GYSUM(F)*DF OVER THE 
APPROPRIATE RANGE OF FREQUENCIES. WRITE OUT THE VARIANCES. 


SUM = 0. 
SUMM = 0. 
DO 900 I = 1, NF 
PRODX = GXSUM(1)*DF 
PRODY = GYSUM(1I)*DF 
SUM = SUM + PRODX 
SUMM = SUMM + PRODY 
900 CONTINUE 
WRITE(FILENO,910) SUM 
WRITE(FILENO,920) SUMM 
910 FORMAT(//,4X,' VARIANCE OF LOD SERIES 
920 FORMAT(//,4X,'VARIANCE OF ATM SERIES 


COMPUTE THE CO-VARIANCE BY INTEGRATING 
WRITE OUT THE CO-VARIANCE. 
5 = 0. 
DO 2000 I = 1, NF 
PRODXY = CXY(1)*DF 
S = St PRODXY 


2000 CONTINUE 


WRITECFILENO,2010) S 


2010 FORMAT(//,4X, ‘COVARIANCE = ',F9.3) 


Sor 
END 


SUBROUTINE TAPER(NTS, X) 


a pe?) 3) 
= ' F9.3) 


imewoO- OPEC LRUMeCKY< EF ). 


MULTIPLYING INPUT DATA BY TAPER FUNCTION 


DIMENSION X(NTS) 
PI = 3. 1415927 
XNTS = NTS 
MA = 0. 1*XNTS 
MB = XNTS-MA 
DO 204 I = 1,NTS 
IF(I. LE.MA) GO TO 201 
IF(1I.GE.MB) GO TO 202 
GO TO 204 
202 XI=NTS-I 
G@ieto 203 
201 XI=I-1 
203 ARG = 5. 0*PI*XI/XNTS 
X(1) = X(1I)*(1. -COS(ARG)*COS( ARG) ) 
204 CONTINUE 
RETURN 
END 


SUBROUTINE MINMAX (A,N, AMIN, AMAX) 
DIMENSION A(1) 


AMIN 
AMAX 


A(1) 
AC 1) 


DO 100 I = 1, N 
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ACT) 


IF (ACI) .LT. AMIN) AMIN 
A(T) 


IF (ACI) .GT. AMAX) AMAX 
100 CONTINUE 
C 


RETURN 
END 
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Fig. 17. Spectral response of 23-point sin? weighted filter: Lowpass filter designed 
primarily to remove solid body tides at 7 and 13.5 days from dLOD¢sxo 
data. 
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Fig. 18. Spectral response of a twice run 75-point sin? filter; Lowpass design per- 
mitted removal of most subseasonal components. Note the abscence of 
side lobes, as well as the relative steepness achieved by running the filter 
twice. 
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Fig. 19. Response of 183-point unweighted moving average filter: Used to isolate 
annual terms. A weighted filter would have required too manv-points, and 
would not have sufficient vertical steepness to include annual frequencies. 
Still, almost 35% of the annual signal is attenuated in order to filter out 


the semiannual terms. 
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Fie, 20. ‘Filter to isolate semiannual terms: Combination of twice run 75-point 
sin? weighted filter and 183-point unweighted subtraction moving average. 
Maximum response is at about 200 davs. 
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Fig. 21. Subseasonal isolation filter; Combination of twice run 75-point weighted 
moving average (see Fig. 19) and 23-point weighted moving average. Peak 
response is at 60 days. 
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